Semesterarbeit - Statistische Datenanalyse DS-C-SD001, ZH-Sa-1, FS25
Sandra Steiner
15.06.2025
library(tidyverse)
library(ggplot2)
library(plotly)
library(GGally)
library(Hmisc)
library(gridExtra)
library(grid)
library(car)
library(viridis)
library(RColorBrewer)
library(dendextend)
library(psych)1 Explorative Datenanalyse
1.1 Laden der Daten
Ich habe die Daten bei World Bank Open Data gefunden und dort heruntergeladen (vgl. Quellenverzeichnis “World Bank Open Data” (n.d.)). Dort habe ich alle Indikatoren für alle Länder für alle erhältlichen Zeiträume abgefragt. Der gesamte Download ist im Ordner “WDI_CSV_2025_01_28” abgelegt (mit Metadaten, i.e. Source, Footnote etc.), welcher bei mir erhältlich gemacht werden kann.
Indicators_comp <- read_csv("WDI_CSV_2025_01_28/WDICSV.csv")head(Indicators_comp)## # A tibble: 6 × 68
## `Country Name` `Country Code` `Indicator Name` `Indicator Code` `1960` `1961`
## <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 Africa Eastern… AFE Access to clean… EG.CFT.ACCS.ZS NA NA
## 2 Africa Eastern… AFE Access to clean… EG.CFT.ACCS.RU.… NA NA
## 3 Africa Eastern… AFE Access to clean… EG.CFT.ACCS.UR.… NA NA
## 4 Africa Eastern… AFE Access to elect… EG.ELC.ACCS.ZS NA NA
## 5 Africa Eastern… AFE Access to elect… EG.ELC.ACCS.RU.… NA NA
## 6 Africa Eastern… AFE Access to elect… EG.ELC.ACCS.UR.… NA NA
## # ℹ 62 more variables: `1962` <dbl>, `1963` <dbl>, `1964` <dbl>, `1965` <dbl>,
## # `1966` <dbl>, `1967` <dbl>, `1968` <dbl>, `1969` <dbl>, `1970` <dbl>,
## # `1971` <dbl>, `1972` <dbl>, `1973` <dbl>, `1974` <dbl>, `1975` <dbl>,
## # `1976` <dbl>, `1977` <dbl>, `1978` <dbl>, `1979` <dbl>, `1980` <dbl>,
## # `1981` <dbl>, `1982` <dbl>, `1983` <dbl>, `1984` <dbl>, `1985` <dbl>,
## # `1986` <dbl>, `1987` <dbl>, `1988` <dbl>, `1989` <dbl>, `1990` <dbl>,
## # `1991` <dbl>, `1992` <dbl>, `1993` <dbl>, `1994` <dbl>, `1995` <dbl>, …Wenn wir das File visuell anschauen, d.h. es entweder in RStudio oder separat als csv öffnen, sehen wir, wie es aufgebaut ist. Wenn wir die verfügbaren Daten für das Jahr 2023 visuell inspizieren und mit den Jahren 2022 und 2021 vergleichen, sehen wir, dass für 2023 die Daten weniger komplett sind. Wir wählen deswegen die Daten für das Jahr 2021, da diese visuell am komplettesten erscheinen.
Indicators_wide_comp <- Indicators_comp %>% select("Country Name", "Indicator Name", "2021") %>% pivot_wider(names_from = "Indicator Name", values_from = "2021")Wir wählen also die Daten fürs 2021 und pivotieren, s.d. die Indikatoren als Spalten angezeigt werden. Nun können wir die Variablen/Indikatoren inspizieren. Ich habe dies schon im Vornherein auf der Website gemacht und ich interessiere mich u.a. für die relative Anzahl Parlamentssitze, welche von Frauen besetzt werden.
1.2 Fragestellungen
Wir können uns folgende Fragen stellen, bevor wir mit der Auswahl der Indikatoren weitermachen.
- Können wir andere Indikatoren finden, von welchen wir annehmen, dass sie mit Gleichberechtigung im Zusammenhang stehen?
- Sind diese Indikatoren mit der relativen Anzahl von Frauen besetzten Parlamentssitzen (linear) korreliert?
- Wenn man solche Indikatoren findet, kann man ein Modell bilden?
- Unterscheiden sich die verschiedenen Länder betreffend dieser von uns gewählten Indikatoren?
- Wenn ja, wie unterscheiden sie sich?
- Kann eine Projektion auf die zwei ersten Hauptkomponenten (PC) eine Antwort auf diese Frage geben?
- Oder können uns verschiedene Clustering-Methoden, Antworten auf die Frage, ob und wie sich die Länder betr. dieser Indikatoren unterscheiden, geben?
Wir wählen Indikatoren aus und überprüfen sie.
df_ind_comp <- Indicators_wide_comp %>% select("Country Name", "GDP per capita (current US$)" , "Female share of employment in senior and middle management (%)" , "Part time employment, female (% of total female employment)" , "Proportion of seats held by women in national parliaments (%)" , "School enrollment, secondary (% gross)","Employment to population ratio, 15+, total (%) (national estimate)", "Government expenditure on education, total (% of GDP)")str(df_ind_comp)## tibble [266 × 8] (S3: tbl_df/tbl/data.frame)
## $ Country Name : chr [1:266] "Africa Eastern and Southern" "Africa Western and Central" "Arab World" "Caribbean small states" ...
## $ GDP per capita (current US$) : num [1:266] 1523 1748 6446 11920 19144 ...
## $ Female share of employment in senior and middle management (%) : num [1:266] NA NA NA NA NA NA NA NA NA NA ...
## $ Part time employment, female (% of total female employment) : num [1:266] NA NA NA NA 25.6 ...
## $ Proportion of seats held by women in national parliaments (%) : num [1:266] 30.5 17.8 18.2 25.4 24.2 ...
## $ School enrollment, secondary (% gross) : num [1:266] 43.8 46.5 69.9 84.4 101.3 ...
## $ Employment to population ratio, 15+, total (%) (national estimate): num [1:266] NA NA NA NA 54.5 ...
## $ Government expenditure on education, total (% of GDP) : num [1:266] 4.77 3.2 NA NA 4.68 ...Um die Daten einfacher handhaben zu können, geben wir den Spalten einfachere Namen.
colnames(df_ind_comp) <- c("Country", "GDP_per_capita", "Females_Management", "Female_parttime", "Females_Parliament","School_enrollment_secondary", "Employment", "Expenditure_education")| Name in Original | neuer Name |
|---|---|
| Country Name | Country |
| GDP per capita (current US$) | GDP_per_capita |
| Female share of employment in senior and middle management (%) | Females_Management |
| Part time employment, female (% of total female employment) | Female_parttime |
| Proportion of seats held by women in national parliaments (%) | Females_Parliament |
| School enrollment, secondary (% gross) | School_enrollment_secondary |
| Employment to population ratio, 15+, total (%) (national estimate) | Employment |
| Government expenditure on education, total (% of GDP) | Expenditure_education |
In dieser Arbeit werden wir den Titel des Indikators bzw. der Variable im Original mit der Bedeutung der Variable gleichsetzen. Für Details zur detaillierten Bedeutung und Erhebung dieser Indikatoren/Variablen verweise ich vollständig auf die Metadaten in WDISeries.csv (können bei Bedarf bei mir erhältlich gemacht werden).
head(df_ind_comp)## # A tibble: 6 × 8
## Country GDP_per_capita Females_Management Female_parttime Females_Parliament
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Africa E… 1523. NA NA 30.5
## 2 Africa W… 1748. NA NA 17.8
## 3 Arab Wor… 6446. NA NA 18.2
## 4 Caribbea… 11920. NA NA 25.4
## 5 Central … 19144. NA 25.6 24.2
## 6 Early-de… 3733. NA 41.5 24.3
## # ℹ 3 more variables: School_enrollment_secondary <dbl>, Employment <dbl>,
## # Expenditure_education <dbl>Wenn wir uns das Resultat anschauen, sehen wir, dass zuoberst verschiedene Regionen/Gruppierungen von Ländern und nicht die einzelnen Länder aufgeführt sind. Wir löschen diesen Part, bzw. speichern die Länder in einem anderen Dataframe
df_ind_country <- df_ind_comp %>% slice(50:266)Wir haben nun ein Dataframe mit 217 Ländern, welche die von uns gewählten Indikatoren aufweisen.
summary(df_ind_country)## Country GDP_per_capita Females_Management Female_parttime
## Length:217 Min. : 214.1 Min. : 7.828 Min. : 6.62
## Class :character 1st Qu.: 2518.6 1st Qu.:27.000 1st Qu.:22.18
## Mode :character Median : 7176.7 Median :34.156 Median :34.92
## Mean : 20303.6 Mean :32.823 Mean :34.29
## 3rd Qu.: 27187.0 3rd Qu.:38.813 3rd Qu.:46.19
## Max. :223823.4 Max. :58.386 Max. :77.84
## NA's :7 NA's :128 NA's :117
## Females_Parliament School_enrollment_secondary Employment
## Min. : 0.00 Min. : 5.46 Min. :27.52
## 1st Qu.:15.06 1st Qu.: 82.87 1st Qu.:51.08
## Median :23.75 Median : 95.23 Median :57.00
## Mean :24.68 Mean : 90.71 Mean :56.02
## 3rd Qu.:33.81 3rd Qu.:104.52 3rd Qu.:61.85
## Max. :61.25 Max. :143.36 Max. :87.28
## NA's :28 NA's :68 NA's :100
## Expenditure_education
## Min. : 0.3819
## 1st Qu.: 3.1151
## Median : 4.2549
## Mean : 4.4708
## 3rd Qu.: 5.4549
## Max. :14.1952
## NA's :52Wir sehen aber, dass es relativ viele NAs hat. Wenn wir alle NAs löschen, erhalten wir ein Dataframe mit 61 Ländern.
df_ind_country_nona <-df_ind_country %>% drop_na()summary(df_ind_country_nona)## Country GDP_per_capita Females_Management Female_parttime
## Length:61 Min. : 828.8 Min. : 7.828 Min. : 6.62
## Class :character 1st Qu.: 4912.6 1st Qu.:28.556 1st Qu.:27.91
## Mode :character Median : 14986.8 Median :34.446 Median :36.18
## Mean : 27212.1 Mean :33.047 Mean :36.74
## 3rd Qu.: 43725.1 3rd Qu.:38.025 3rd Qu.:49.10
## Max. :133711.8 Max. :57.244 Max. :77.84
## Females_Parliament School_enrollment_secondary Employment
## Min. : 0.00 Min. : 41.86 Min. :31.00
## 1st Qu.:20.80 1st Qu.: 91.50 1st Qu.:51.08
## Median :28.32 Median :101.83 Median :56.35
## Mean :30.00 Mean :101.15 Mean :55.93
## 3rd Qu.:40.00 3rd Qu.:111.80 3rd Qu.:59.40
## Max. :61.25 Max. :143.36 Max. :79.19
## Expenditure_education
## Min. :1.507
## 1st Qu.:3.747
## Median :4.730
## Mean :4.762
## 3rd Qu.:5.497
## Max. :8.218str(df_ind_country_nona)## tibble [61 × 8] (S3: tbl_df/tbl/data.frame)
## $ Country : chr [1:61] "Albania" "Angola" "Argentina" "Austria" ...
## $ GDP_per_capita : num [1:61] 6413 1926 10738 53649 7490 ...
## $ Females_Management : num [1:61] 30.7 11 37.4 34.3 44 ...
## $ Female_parttime : num [1:61] 26.9 34.8 57.5 65 22.2 ...
## $ Females_Parliament : num [1:61] 35.7 29.5 44.7 40.4 40 ...
## $ School_enrollment_secondary: num [1:61] 97.7 52.5 116.6 102.3 95.2 ...
## $ Employment : num [1:61] 52.8 64.3 55.1 57.2 67.3 ...
## $ Expenditure_education : num [1:61] 3.02 2.3 4.64 5.49 4.61 ...Sehen wir uns zuerst einmal die einzelnen Variablen an. Da viele der im Folgenden angewendeten Methoden einen vollständigen Datensatz erfordern, werden wir im Anschluss jeweils mit den Ländern weiterarbeiten, welche für alle gewählten Indikatoren Einträge aufweisen.
1.3 Analyse der einzelnen Variablen
str(df_ind_comp)## tibble [266 × 8] (S3: tbl_df/tbl/data.frame)
## $ Country : chr [1:266] "Africa Eastern and Southern" "Africa Western and Central" "Arab World" "Caribbean small states" ...
## $ GDP_per_capita : num [1:266] 1523 1748 6446 11920 19144 ...
## $ Females_Management : num [1:266] NA NA NA NA NA NA NA NA NA NA ...
## $ Female_parttime : num [1:266] NA NA NA NA 25.6 ...
## $ Females_Parliament : num [1:266] 30.5 17.8 18.2 25.4 24.2 ...
## $ School_enrollment_secondary: num [1:266] 43.8 46.5 69.9 84.4 101.3 ...
## $ Employment : num [1:266] NA NA NA NA 54.5 ...
## $ Expenditure_education : num [1:266] 4.77 3.2 NA NA 4.68 ...Für jedes Land haben wir das Bruttoinlandprodukt (BIP) pro Kopf (GDP per Capita in US$) (“CDP_per_capita”), den Anteil weiblicher Personen bei Anstellungen im höheren und mittleren Management (“Females_Management), den Anteil weiblicher Teilzeitarbeit an der gesamten weiblichen Arbeit (Female_parttime), den Anteil an weiblichen Parlamentsmitgliedern, den Prozentsatz von Einschulungen in die zweite Schulstufe (School_enrollment_secondary), Verhältnis von Anstellungen zur Bevölkerungszahl ab 15 Jahre (Employment) und die Staatsausgaben für die Bildung als Prozentsatz des BIPs (Expenditure_education).
Für alle diese Variablen erstellen wir jeweils 4 Plots zwecks explorativer Datenanalyse und interpretieren diese.
var<-"GDP_per_capita"
mu <- mean(df_ind_country[[var]], na.rm = TRUE)
sigma <- sd(df_ind_country[[var]], na.rm = TRUE)
plot1 <-ggplot(data=df_ind_country)+
geom_point(mapping=aes_string(x="Country",y=var))+
theme_bw()+
theme(axis.text.x = element_blank())
plot2 <- ggplot(data=df_ind_country)+
geom_histogram(mapping=aes_string(x=var), binwidth=2500)+
theme_bw()
plot3 <- ggplot(data=df_ind_country)+
geom_boxplot(mapping=aes_string(y=var),fill="lightblue")+
theme_bw()+
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank())
plot4 <- ggplot(data=df_ind_country, aes(sample = !!sym(var))) +
geom_qq() +
geom_qq_line(color="red")+
theme_bw()
list_of_plots <- list(plot1,plot2, plot3, plot4)
text <- paste("Variable =", var)
grid.arrange(arrangeGrob(grobs=list_of_plots, ncol=2), textGrob(text, gp = gpar(fontsize = 12, fontface ="bold")), heights = c(1, 0.2))
Sehr viele Länder haben offenbar ein relativ niedriges BIP pro Kopf und einige wenige ein sehr hohes. Dieser Indikator ist eindeutig nicht normalverteilt. Der Boxplot und der QQ-Plot unterstreichen dieses Bild. Die beiden Enden des QQ-Plots sind nach oben gebogen im Vergleich zur Geraden, welche bei einer Normalverteilung entstehen würde.
var<-"Females_Management"
plot1 <-ggplot(data=df_ind_country)+
geom_point(mapping=aes_string(x="Country",y=var))+
ylim(0,100)+
theme_bw()+
theme(axis.text.x = element_blank())
plot2 <- ggplot(data=df_ind_country)+
geom_histogram(mapping=aes_string(x=var))+
xlim(0,100)+
theme_bw()
plot3 <- ggplot(data=df_ind_country)+
geom_boxplot(mapping=aes_string(y=var),fill="lightblue")+
ylim(0,100)+
theme_bw()+
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank())
plot4 <- ggplot(data=df_ind_country, aes(sample = !!sym(var))) +
geom_qq() +
geom_qq_line(color="red")+
theme_bw()
list_of_plots <- list(plot1,plot2, plot3, plot4)
text <- paste("Variable =", var)
grid.arrange(arrangeGrob(grobs=list_of_plots, ncol=2), textGrob(text, gp = gpar(fontsize = 12, fontface ="bold")), heights = c(1, 0.2))
Der Prozentsatz der weiblichen Angestellten im mittleren und höheren Management scheint nahezu normalverteilt zu sein. Der Median und auch der Mittelwert scheinen dabei ungefähr bei 35 Prozent zu liegen. Man sieht auch von Auge, dass die Punkte auf dem Scatterplot schön gleichmässig um den Mittelwert streuen. Auf dem QQ-Plot sieht man, dass die Normalverteilung gut angenähert wird, nur an den Enden, wo aber nur wenige Punkte sind, gibt es Abweichungen von der Referenzgeraden. Die Enden sind ein evt. ein bisschen “schwer”, es sind häufiger extremere Werte als bei einer Normalverteilung möglich.
var<-"Female_parttime"
plot1 <-ggplot(data=df_ind_country)+
geom_point(mapping=aes_string(x="Country",y=var))+
ylim(0,100)+
theme_bw()+
theme(axis.text.x = element_blank())
plot2 <- ggplot(data=df_ind_country)+
geom_histogram(mapping=aes_string(x=var))+
xlim(0,100)+
theme_bw()
plot3 <- ggplot(data=df_ind_country)+
geom_boxplot(mapping=aes_string(y=var),fill="lightblue")+
ylim(0,100)+
theme_bw()+
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank())
plot4 <- ggplot(data=df_ind_country, aes(sample = !!sym(var))) +
geom_qq() +
geom_qq_line(color="red")+
theme_bw()
list_of_plots <- list(plot1,plot2, plot3, plot4)
text <- paste("Variable =", var)
grid.arrange(arrangeGrob(grobs=list_of_plots, ncol=2), textGrob(text, gp = gpar(fontsize = 12, fontface ="bold")), heights = c(1, 0.2))
Der Anteil von Teilzeitangestellten an allen weiblichen Angestellten ist ebenfalls um einen Median von ca. 35 Prozent verteilt. Es ist aber insbesondere im QQ-Plot eine gewisse Linksschiefe zu sehen. Es ist am linken Rand eine deutliche Abweichung von der Geraden nach oben zu sehen. Die Verteilung hat weniger extreme, kleine Werte als die Normalverteilung.
var<-"Females_Parliament"
plot1 <-ggplot(data=df_ind_country)+
geom_point(mapping=aes_string(x="Country",y=var))+
ylim(0,100)+
theme_bw()+
theme(axis.text.x = element_blank())
plot2 <- ggplot(data=df_ind_country)+
geom_histogram(mapping=aes_string(x=var))+
xlim(0,100)+
theme_bw()
plot3 <- ggplot(data=df_ind_country)+
geom_boxplot(mapping=aes_string(y=var),fill="lightblue")+
ylim(0,100)+
theme_bw()+
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank())
plot4 <- ggplot(data=df_ind_country, aes(sample = !!sym(var))) +
geom_qq() +
geom_qq_line(color="red")+
theme_bw()
list_of_plots <- list(plot1,plot2, plot3, plot4)
text <- paste("Variable =", var)
grid.arrange(arrangeGrob(grobs=list_of_plots, ncol=2), textGrob(text, gp = gpar(fontsize = 12, fontface ="bold")), heights = c(1, 0.2))
Der Anteil von weiblichen Mitgliedern in den Parlamenten der erfassten Ländern weist einen Median von ca. 25 Prozent auf. Auch hier sind die Datenpunkte relativ gleichmässig um den Median verteilt. Auf dem QQ-Plot sieht man, dass die Gerade ausser am äussersten linken Rand gut angenähert wird. Am linken Rand sieht man, dass extreme Werte weniger häufig als bei der Normalverteilung sind.
var<-"School_enrollment_secondary"
plot1 <-ggplot(data=df_ind_country)+
geom_point(mapping=aes_string(x="Country",y=var))+
ylim(0,100)+
theme_bw()+
theme(axis.text.x = element_blank())
plot2 <- ggplot(data=df_ind_country)+
geom_histogram(mapping=aes_string(x=var))+
xlim(0,100)+
theme_bw()
plot3 <- ggplot(data=df_ind_country)+
geom_boxplot(mapping=aes_string(y=var),fill="lightblue")+
ylim(0,100)+
theme_bw()+
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank())
plot4 <- ggplot(data=df_ind_country, aes(sample = !!sym(var))) +
geom_qq() +
geom_qq_line(color="red")+
theme_bw()
list_of_plots <- list(plot1,plot2, plot3, plot4)
text <- paste("Variable =", var)
grid.arrange(arrangeGrob(grobs=list_of_plots, ncol=2), textGrob(text, gp = gpar(fontsize = 12, fontface ="bold")), heights = c(1, 0.2))
Der Anteil der Personen, welche in eine Schule der sekundären Bildungsstufe eingeschrieben waren, ist stark rechtsschief verteilt. Auf dem Scatterplot sieht man, dass die Punkte nicht gleichmässig um einen Mittelwert streuen. Der QQ-Plot zeigt ebenfalls, dass die Daten stark von einer Normalverteilung abweichen.
var<-"Employment"
plot1 <-ggplot(data=df_ind_country)+
geom_point(mapping=aes_string(x="Country",y=var))+
ylim(0,100)+
theme_bw()+
theme(axis.text.x = element_blank())
plot2 <- ggplot(data=df_ind_country)+
geom_histogram(mapping=aes_string(x=var))+
xlim(0,100)+
theme_bw()
plot3 <- ggplot(data=df_ind_country)+
geom_boxplot(mapping=aes_string(y=var),fill="lightblue")+
ylim(0,100)+
theme_bw()+
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank())
plot4 <- ggplot(data=df_ind_country, aes(sample = !!sym(var))) +
geom_qq() +
geom_qq_line(color="red")+
theme_bw()
list_of_plots <- list(plot1,plot2, plot3, plot4)
text <- paste("Variable =", var)
grid.arrange(arrangeGrob(grobs=list_of_plots, ncol=2), textGrob(text, gp = gpar(fontsize = 12, fontface ="bold")), heights = c(1, 0.2))
Der Anteil Beschäftigter bei über 15-jährigen Personen streut um einen Median von ca. 60 Prozent. Die Verteilung ähnelt zwar einer Normalverteilung, aber es scheint, dass die Spitze extremer zuläuft und dafür an den äusseren Rändern extreme Werte häufiger sind als bei einer Normalverteilung (schwere Enden).
var<-"Expenditure_education"
plot1 <-ggplot(data=df_ind_country)+
geom_point(mapping=aes_string(x="Country",y=var))+
ylim(0,10)+
theme_bw()+
theme(axis.text.x = element_blank())
plot2 <- ggplot(data=df_ind_country)+
geom_histogram(mapping=aes_string(x=var), bins=30)+
xlim(0,10)+
theme_bw()
plot3 <- ggplot(data=df_ind_country)+
geom_boxplot(mapping=aes_string(y=var),fill="lightblue")+
ylim(0,10)+
theme_bw()+
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank())
plot4 <- ggplot(data=df_ind_country, aes(sample = !!sym(var))) +
geom_qq() +
geom_qq_line(color="red")+
theme_bw()
list_of_plots <- list(plot1,plot2, plot3, plot4)
text <- paste("Variable =", var)
grid.arrange(arrangeGrob(grobs=list_of_plots, ncol=2), textGrob(text, gp = gpar(fontsize = 12, fontface ="bold")), heights = c(1, 0.2))
Der Anteil der Staatsausgaben für Bildung gemessen am BIP ist für die Länder um einen Median von ca. 4 Prozent verteilt. Hier musste die x-Achse, welche für die anderen Plots auf 0 bis 100 Prozent fixiert war, auf 0 bis 10 Prozent eingeschränkt werden, um eine sinnvolle Ansicht zu erhalten. Auch hier liegt keine perfekte Normalverteilung vor, der QQ-Plot zeigt Abweichungen nach oben am rechten Ende, was bedeutet, dass dort extremere Werte häufiger als bei einer Normalverteilung vorkommen.
1.3.1 Fazit
Das BIP (GDP) ist deutlich nicht normalverteilt. Gleiches gilt für die Variable “School_enrollment_secondary”. Ansonsten sind die Prozentsätze jeweils näherungsweise normalverteilt. Dies sollte uns nicht erstaunen, wissen wir doch, dass wenn die Stichprobe genug gross und der Anteilswert nicht zu extrem ist, die Normalverteilung die Binomialverteilung approximiert (vgl. (Diez, Çetinkaya-Rundel, and Barr 2019)).
1.4 Korrelationen
Sehen wir uns die Korrelationen und Scatterplots zwischen den Variablen an. Hier wird der Pearson-Korrelationskoeffizient gezeigt.
# Streudiagramm-Matrix mit Korrelationen
ggpairs(
df_ind_country[, colnames(df_ind_country_nona[, -1])],
title = "Streudiagramm-Matrix mit Pearson-Korrelationskoeffizient",
upper = list(continuous = "cor"),
lower = list(continuous = "points"),
diag = list(continuous = "densityDiag")
)
Wenn wir uns an die Fragestellungen erinnern, dann sehen wir, dass vor allem mit dem GDP per Capita zwar ein Zusammenhang mit anderen Variablen existiert, aber nicht unbedingt ein linearer. Wenn wir uns die Variable Frauenanteil im Parlament anschauen, gibt es dort eher (schwach) lineare Zusammenhänge mit anderen Variablen. Insbesondere ist der Korrelations-Koeffizient zwischen Females_parliament und Female_parttime 0.405.
Wenn wir später ein Modell bilden möchten, müssen wir für jedes Land alle uns interessierenden Daten haben. Schauen wir uns deswegen den Plot noch einmal an, aber nur mit den 61 Ländern, für welche alle Daten vorhanden sind.
# Streudiagramm-Matrix mit Korrelationen
ggpairs(
df_ind_country_nona[, colnames(df_ind_country_nona[, -1])],
title = "Streudiagramm-Matrix mit Pearson-Korrelationskoeffizient",
upper = list(continuous = "cor"),
lower = list(continuous = "points"),
diag = list(continuous = "densityDiag")
)
Für diese Länder ist die Korrelation zwischen Females_parliament und Female_parttime 0.515 und zwischen Female_parliament und Ausgaben für Bildung 0.497.
Diese zwei Variablen werden wir später für die Modellbildung im Kapitel 2 auswählen.
Selektieren wir den Datensatz also nach diesen Variablen und schliessen die NAs aus.
df_ind_modell <-df_ind_country %>% select(Females_Parliament, Female_parttime, Expenditure_education) %>% drop_na()Dieser Datensatz beinhaltet nun 86 Länder. Schauen wir uns noch einmal die Korrelationen an:
# Streudiagramm-Matrix mit Korrelationen
ggpairs(
df_ind_modell[, ],
title = "Streudiagramm-Matrix mit Pearson-Korrelationskoeffizient",
upper = list(continuous = "cor"),
lower = list(continuous = "points"),
diag = list(continuous = "densityDiag")
)
Wir sehen, dass es die Variablen nun eine kleinere Korrelation aufweisen. Da wir aber keine besseren Daten haben, schauen wir die Modellbildung mit diesen Daten an. Die Korrelation zwischen den beiden Variablen Female_parttime und Expenditure_education ist 0.3. Diesen Wert werden wir später bei der Prüfung der Multikollinerarität noch einmal anschauen.
2 Multiple lineare Regression
2.1 Multiple lineare Regression - Theoretische Grundlagen
Bei der multiplen linearen Regression geht man davon aus, dass die Abhängigkeit einer Variablen \(Y\) von erklärenden Variablen \(X_1,...,X_p\) linear ist.
Das heisst es gibt Parameter \({\beta}_0,...,{\beta}_p\), s.d.:
\[Y={\beta}_0+{\beta}_1X_1+...+{\beta}_pX_p + \epsilon \]
Die Störgrösse \(\epsilon\) fasst dabei den möglichen Einfluss von weiteren unbekannten Variablen zusammen und wir gehen davon aus, dass dieser Einfluss additiv wirkt (Handl and Kuhlenkasper 2017).
Ist \(p=1\) handelt spricht man von linearer Einfachregression, ansonsten von der multiplen linearen Regression (Handl and Kuhlenkasper 2017).
Da die Parameter \({\beta}_0,...,{\beta}_p\) nicht im Vornherein bekannt sind, müssen sie aus den Daten geschätzt werden. Man macht dies mit der Kleinsten-Quadrate-Methode (Handl and Kuhlenkasper 2017).
Man hat in den Daten also Variablen \(Y\) und \(X_1,...,X_p\) mit ihren jeweiligen Realisierungen zu Messzeitpunkten \(1\),…,\(n\), welche wir mit \(y_i\), \(x_{ik}\), für \(i=1,...,n\); \(k=1,...,p\) bezeichnen.
Wir werden der Einfachheit halber noch \(x_{0k}=1\) einführen, für alle \(k=1,...,p\). Wir suchen dann die \({\beta}_0,...,{\beta}_p\), s.d.
\[ \sum_{i=1}^n (y_i - \hat{y}_i)^2 = \sum_{k=1}^n \left(y_i - \sum_{k=0}^p \beta_k x_{ik} \right)^2 \]
so klein wie möglich ist. Wir legen also eine Hyperebene so durch unseren Punkteraum, dass die Summe der quadratischen Abstände zur Ebene minimal ist (Handl and Kuhlenkasper 2017).
Ebenso unbekannt ist der Einfluss der Störgrösse. Hat man die \({\beta}_0,...,{\beta}_p\) so festgelegt, dass die Summe der quadratischen Abstände zur Ebene minimal ist, kann man die Residuen berechnen. Also für jedes \(y_i\) berechnen wir den Abstand zu \(\hat{y}_i={\beta}_0+{\beta}_1x_1+...+{\beta}_px_p\), indem wir \(e_i=y_i-\hat{y}_i\) berechnen. Die Abstände \(e_i\) nennt man Residuen.
Die Varianz \(\sigma^2\) der Störgrösse \(\epsilon\), mit \(\epsilon = (\epsilon_1, ..., \epsilon_n)\) wird mittels der Residuen \(e_i\) geschätzt. Als Schätzer wird eine leicht abgewandelte Version der Stichprobenvarianz der Residuen verwendet:
\[\sigma^2=\frac{1}{n-p-1}\sum_1^n e_i^2\]
Für die Herleitung verweise ich vollständig auf den entsprechenden Abschnitt im Buch von Handl und Kuhlenkasper (Handl and Kuhlenkasper 2017).
Damit unsere Daten in sinnvoller Weise mit multipler linearer Regression modelliert werden können, müssen verschiedene Voraussetzungen erfüllt sein. In der Regel nehmen wir an, dass sie diese erfüllt sind und überprüfen dies im Detail nach der Modellbildung.
Es handelt sich dabei um folgende Voraussetzungen (zitiert aus dem Tutorial von Markus Geuss der Schulungsunterlagen):
- Lineare Beziehung zwischen den Prädiktoren und der abhängigen Variablen
- Homoskedastizität (konstante Varianz der Fehler/Residuen)
- Unabhängigkeit der Fehler/Residuen (keine Autokorrelation)
- Normalverteilung der Residuen
- Keine Multikollinearität der Prädiktoren/erklärenden Variablen
Abhängig davon, wie gut diese Voraussetzungen erfüllt sind, lässt sich dann auch die Modellgüte beurteilen.
Dafür zieht man bei der multiplen linearen Regression zusätzlich zum Bestimmtheitsmass \(R^2\)
\[ R^2 = 1 - \frac{SS_{\text{residual}}}{SS_{\text{total}}} = \frac{\sum_{i} \left( \hat{y}_i - \overline{y} \right)^2}{\sum_{i} \left( y_i - \overline{y} \right)^2} \]
auch adjusted \(R^2\) hinzu:
\[ R^2_{\text{adjusted}} = 1 - \left( \frac{SS_{\text{residual}} / (n - p - 1)}{SS_{\text{total}} / (n - 1)} \right) \] \(R^2\) gibt dabei das Verhältnis von erklärter zu gesamter Varianz wieder. Bei adjusted \(R^2\) handelt es sich um eine um die Anzahl der Variablen korrigierte Version, welche im Falle von mehreren Variablen aussagekräftiger ist (vgl.Tutorial von Markus Geuss der Schulungsunterlagen).
Die Voraussetzungen und die Beurteilung der Modellgüte werden wir im Folgenden direkt anhand der Daten weiterbesprechen.
2.2 Multiple lineare Regression - Modellbildung
Wir erstellen ein Modell für eine multiple lineare Regression um den Zusammenhang zwischen Anzahl weiblicher Parlamentsmitglieder und den Ausgaben für die Bildung, sowie dem Anteil weiblicher Teilzeitangestellter an den weiblichen Angestellten im Land zu modellieren. Wir nehmen als Grundlage den Datensatz mit 86 Ländern, welcher für alle vorkommenden Länder die relevanten Variablen aufweist.
# Erstellung des Modells
modell <- lm(Females_Parliament ~ Female_parttime + Expenditure_education, data = df_ind_modell)
# Zusammenfassung des Modells anzeigen
summary(modell)##
## Call:
## lm(formula = Females_Parliament ~ Female_parttime + Expenditure_education,
## data = df_ind_modell)
##
## Residuals:
## Min 1Q Median 3Q Max
## -30.4526 -6.0303 -0.8618 6.4128 25.6913
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.77609 4.21140 1.846 0.06839 .
## Female_parttime 0.26875 0.07659 3.509 0.00073 ***
## Expenditure_education 2.43917 0.84823 2.876 0.00512 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.5 on 83 degrees of freedom
## Multiple R-squared: 0.2595, Adjusted R-squared: 0.2416
## F-statistic: 14.54 on 2 and 83 DF, p-value: 3.856e-062.2.1 Prüfen der Modellvoraussetzungen
Wir können nun die Modellvoraussetzungen prüfen:
- Lineare Beziehung zwischen den Prädiktoren und der abhängigen Variablen
Wir haben in den Scatterplots mit den Korrelationen bereits gesehen, dass der Pearson-Korrelationskoeffizient mit ca. 0.39 und ca. 0.43 eher tief ist. Wenn wir uns die Komponent-Residuen-Plots ansehen, müssten die Datenpunkte bei perfekter Linearität um die blaue Gerade streuen. Die Plots zeigen für eine erklärende Variable (Prädiktor) \(X_i\) die Komponente \(\beta_ix_{ik}\) plus das Residual \(e_k\) gegen die jeweiligen \(x_{ik}\). Das heisst, die anderen Prädiktoren werden sozusagen heraus gerechnet und wir sehen nur den Einfluss des Prädiktors \(X_i\), welcher linear sein müsste.
crPlots(modell)
Wir sehen, was wir schon bei der Korrelationsmatrix beobachtet haben, der lineare Zusammenhang ist nicht gegeben und die Punkte streuen in einer Wolke um die Gerade.
- Homoskedastizität (konstante Varianz der Fehler/Residuen)
Auch hier können wir die diagnostischen Plots anschauen. Wenn die Varianz konstant wäre, würden die standardisierten Residuen gleichmässig um eine horizontale Gerade streuen. Wir sehen, dass es hier eine Abweichung davon gibt, i.e. es gibt einen Knick in den Plots der Residuen. Während die Varianz am Anfang zwischen 0 und 30 (Fitted values) eher grösser ist, also die standardisierten Residuen mehr streuen, streuen sie nach 30 (Fitted values) eher weniger. Ausserdem scheinen die (standardisierten) Residuen nach 30 auch kleiner zu werden. Die Varianz der Residuen ist also eher nicht konstant.
plot(modell, which = 3)
Wir können uns dazu auch noch einen statistischen Test ansehen. Dieser hat als Nullhypothese, dass die Varianz konstant ist.
ncvTest(modell)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.09485084, Df = 1, p = 0.7581Hier kann die Nullhypothese aufgrund des Tests zwar nicht verworfen werden, i.e. p ist nicht kleiner als 0.05, aber wir haben vorher bereits gesehen, dass es hier ein Problem gibt und eher keine Homoskedastizität vorliegt.
- Unabhängigkeit der Fehler/Residuen (keine Autokorrelation)
Es wird auch vorausgesetzt bzw. angenommen, dass die Residuen unabhängig voneinander sind. Wenn dies nicht der Fall ist, kann das zu verzerrten Schätzungen führen. Man spricht von Autokorrelation, da der Fehler eines Datenpunktes mit den Fehlern der anderen Datenpunkte korreliert.
durbinWatsonTest(modell)## lag Autocorrelation D-W Statistic p-value
## 1 -0.01400003 2.005845 0.994
## Alternative hypothesis: rho != 0Der Durbin-Watson Wert (D-W Statistic) sollte nahe bei 2 liegen, wenn keine Autokorrelation vorliegt. Der Wert für die Autokorrelation sollte nahe bei Null liegen. Das ist hier der Fall. Der p-Wert ist hoch, was ebenfalls die Nullhypothese (keine Autokorrelation) stützt. Der Test zeigt also, dass hier die Annahme der Unabhängigkeit der Residuen nicht verletzt ist.
- Normalverteilung der Residuen
Für die Überprüfung dieser Bedingung können wir den QQ-Plot heranziehen.
plot(modell, which = 2)
Der QQ-Plot der standardisierten Residuen sieht grundsätzlich nicht schlecht aus. An den Rändern sieht man leichte Abweichungen, was jedoch auch zu einem gewissen Grad normal ist, da man dort weniger “Sampling” hat. Die Normalverteilung der Residuen scheint nicht verletzt zu sein.
- Keine Multikollinearität der Prädiktoren
Wir haben hier nur zwei Prädiktoren/erklärende Variablen, “Female_parttime” und “Expenditure_education”. Im Abschnitt 1.4. Korrelationen haben wir gesehen, dass die Korrelation (Pearson) zwischen den beiden Variablen ca. 0.3 beträgt. Damit ist diese Voraussetzung in unserem Beispiel eher unproblematisch. Wir können noch den VIF (Variance Inflation Factor) bestimmen, welcher für eine erklärende Variable \(X_i\) folgendermassen berechnet wird:
\[ VIF_i=\frac{1}{1-R_i^2}\] wobei \(R_i^2\) das Bestimmtheitsmass ist, welches errechnet wird, wenn man \(X_i\) mit Hilfe einer multiplen linearen Regression der übrigen erklärenden Variablen modelliert.
vif(modell)## Female_parttime Expenditure_education
## 1.095614 1.095614Wenn der VIF-Wert für eine Variable \(X_i\) nahe bei \(1\) ist, bedeutet das, dass \(R_i\) nahe bei \(0\) ist und somit nicht (gut) linear durch die restlichen Variablen modelliert/erklärt werden kann. Das heisst, es liegt keine problematische Multikollinearität vor.
2.2.2 Güte des Modells
Schauen wir uns das Summary des Modells an.
summary(modell)##
## Call:
## lm(formula = Females_Parliament ~ Female_parttime + Expenditure_education,
## data = df_ind_modell)
##
## Residuals:
## Min 1Q Median 3Q Max
## -30.4526 -6.0303 -0.8618 6.4128 25.6913
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.77609 4.21140 1.846 0.06839 .
## Female_parttime 0.26875 0.07659 3.509 0.00073 ***
## Expenditure_education 2.43917 0.84823 2.876 0.00512 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.5 on 83 degrees of freedom
## Multiple R-squared: 0.2595, Adjusted R-squared: 0.2416
## F-statistic: 14.54 on 2 and 83 DF, p-value: 3.856e-06Die Koeffizienten sind grundsätzlich signifikant, ausser das Intercept, welches das Signifikanzlevel von 0.05 knapp nicht erreicht.
Der F-Wert ist ebenfalls signifikant. Das heisst die Nullhypothese (alle \(\beta_i=0\) ausser \(\beta_0\) (Intercept)) kann verworfen werden und das Modell erklärt signifikant mindestens einen Teil der Varianz.
Das adjusted \(R^2\) ist jedoch relativ niedrig, genauso das \(R^2\). Das Bestimmtheitsmass \(R^2\) ist das Verhältnis zwischen erklärter Varianz und gesamter Varianz, es gibt also den Anteil der Varianz an, welcher durch das Modell erklärt wird. Bei mehreren erklärenden Variablen steigt dies jedoch automatisch pro zugefügter Variable an, ohne dass sicher ist, dass diese Variable wirklich etwas zur Erklärung der Varianz beiträgt. Adjustet \(R^2\) ist deswegen eine um diesen Einfluss von multiplen erklärenden Variablen korrigierte Version. In unserem Fall sehen wir dass lediglich ca. 25 Prozent der Varianz in der Variable Females_Parliament durch unser Modell erklärt werden.
(Bemerkung: Die Plots und Anmerkungen zu den Abschnitten “Prüfen der Modellvoraussetzungen” und “Modellgüte” basieren zu weiten Teilen auf dem Tutorial von Markus Geuss, welches Teil der Unterrichtsmaterialen war.)
2.2.3 Fazit
Wir haben gesehen, dass mindestens zwei der Voraussetzungen nicht erfüllt waren. Insbesondere natürlich die wichtigste, die Linearität. Es war schon nach der explorativen Datenanalyse klar, dass zwischen den zwei erklärenden Variablen und der abhängigen Variablen höchstens ein sehr schwacher linearer Zusammenhang bestand. Zudem haben wir gesehen, dass lediglich ca. 25 Prozent der Varianz in den Daten durch unser Modell erklärt werden.
Jedoch kann je nach Anwendungsfall ein solches Modell trotzdem eingesetzt werden, auch wenn nicht alle Voraussetzungen erfüllt sind und die Modellgüte nicht optimal ist. Dies kann beispielsweise dann sinnvoll sein, wenn man in einem grösseren Datensatz fehlende Werte schätzen und den Datensatz damit ergänzen möchte.
3 PCA
Im Teil 2 haben wir uns mit der Fragestellung, ob wir eine Variable in linearer Abhängigkeit von anderen Variablen modellieren können, befasst. Das nächste Kapitel könnte uns bei diesem Thema ebenfalls helfen. Die PCA (Principal Component Analysis - oder auf Deutsch Hauptkomponentenanalyse) kann bei Multikollinearität in den Daten eine Lösung bieten.
Wir werden hier aber mit Hilfe der PCA Folgendes untersuchen: Unterscheiden sich die verschiedenen Länder betreffend der von uns gewählten Indikatoren? Wenn ja, wie unterscheiden sie sich? Kann eine Projektion auf die zwei ersten Hauptkomponenten (PC) eine Antwort auf diese Frage geben? Bzw. können wir anhand einer Projektion auf die Hauptkomponenten ein Muster erkennen, insbesondere wenn wir die Zugehörigkeit der Länder zu verschiedenen geografischen Weltregionen miteinbeziehen?
3.1 PCA - Theoretische Grundlagen
Oft haben wir einen grossen Datensatz mit vielen Variablen, zum Teil können diese korreliert sein, wie oben bereits erwähnt. Nun gibt es mehrere Ansätze, warum eine Reduktion dieses Datensatzes wünschenswert ist. Einerseits um eine Modellierung zu vereinfachen oder zu verbessern und korrelierte Variablen aus dem Datensatz zu entfernen, ohne Informationen zu verlieren. Andererseits möchten wir vielleicht die einzelnen gemessenen Objekte (hier Länder) “nach Grösse” anordnen oder als Scatterplot darstellen. Dies aber, indem wir alle Variablen einbeziehen (Handl and Kuhlenkasper 2017).
Wir können dies erreichen, indem wir den Datensatz auf die sogenannten Principal Components (Hauptkompenten) reduzieren. Projektion der Daten auf die erste Hauptkomponente liefert dann eine eindimensionale Anordnung der Objekte und Projektion auf die ersten zwei Hauptkomponenten liefert einen Scatterplot, welcher wesentliche Strukturen und Zusammenhänge in den Daten aufzeigen kann.
Die erste Hauptkomponente ist dabei die Richtung, in welcher die Daten am meisten Varianz aufweisen. Die zweite Hauptkomponente ist orthogonal zur ersten und erfasst von allen orthogonalen Richtungen diejenige, welche wiederum am meisten Varianz aufweist. Dieses Prinzip setzt sich so fort, wobei die Varianz, welche die Hauptkomponenten jeweils für sich genommen erklären, stetig abnimmt. Hat es \(p\) Variablen gibt es \(p\) Hauptkomponenten.
Formal lässt sich dies so zusammenfassen:
(Ich zitiere hier im Folgenden gekürzt aus (James et al. 2021) und verwende teilweise auch die Matrixnotation aus (Handl and Kuhlenkasper 2017).)
Wir haben \(n\) Beobachtungen/Objekte und wir haben \(p\) Variablen \(X_1,...,X_p\), welche wir für jedes Objekt gemessen haben. Wir haben also eine Datenmatrix \(X\) aus \(\mathbb{R}^{n\times p}\)
Die erste Hauptkomponente ist die normalisierte Linearkombination der Variablen:
\[ Z_1=a_{11}X_1+a_{21}X_2+ ... +a_{p1}Xp \]
welche die grösste Varianz hat. Normalisiert heisst, dass \(\sum_{j=1}^p{a_{j1}^2=1}\).
Die Koeffizienten \(a_{j1}\) nennen wir die Loadings der ersten Hauptkomponente. Wir können dies auch als Vektor schreiben, \(a_1=(a_{11}, ..., a_{p1})^T\). \(a_1\) gibt die Richtung der Hauptkomponente im alten Variablenraum an.
Für \(X\) können wir annehmen, dass \(X\) um den den Nullpunkt zentriert ist. Dadurch haben alle p Variablen den Mittelwert \(0\).
Es handelt sich um ein Optimierungsproblem, bei welchem
\[ \text{max}_{a_1} \text{Var}(Xa_1) \] unter der Bedingung \(a_1^Ta_1=1\) optimiert wird.
Es stellt sich heraus, dass \(a_1\) ein Eigenvektor der Kovarianzmatrix \[ \Sigma = \frac{1}{n-1} X^\top X \]
ist. Der Eigenwert \(\lambda_1\) zu \(a_1\) ist dabei genau die Varianz in Richtung der ersten Hauptkomponente:
\[ \lambda_1=\text{Var}(Z_1)=\text{Va}r(Xa_1) \]
Für die zweite Hauptkomponente \(Z_2=a_{12}X_1+...+a_{p2}X_p\) muss gelten, dass \(a_2^Ta_2=1\) und \(a_2^Ta_1=0\), also die Richtungen müssen orthonormal zueinander sein. Wiederum muss die verbleibende Varianz maximiert werden.
Dieses Verfahren wird fortgesetzt. Es gilt, dass die \(a_i\) jeweils die Eigenvektoren der der Kovarianzmatrix \(\Sigma\) zu Eigenwerten \(\lambda_i=Var(Xa_i)\) sind. Wobei aufgrund des aufbauenden Verfahrens auch gilt \(\lambda_1 \geq \lambda_2 \geq... \geq \lambda_p\).
Wir erhalten durch Multiplikation unserer Datenmatrix mit den Loadings der Hauptkomponenten, \(A\in \mathbb{R}^{p x p}=[a_1,a_2,...,a_p]\):
\[ XA=Z \] wobei \(Z \in \mathbb{R}^{n\times p}\) die Matrix der neuen Koordinaten unserer Beobachtungen \(1\) bis \(n\) im Hauptkomponentenraum bildet.
Meistens wird dann jeweils auf einige wenige, e.g. \(k<p\) Hauptkomponenten reduziert, welche zusammen schon einen Grossteil der Varianz in den Daten erklären. Man reduziert \(A\) dann auf \(A^k \in \mathbb{R}^{p\times k}\), mit \(A^k=[a_1, ..., a_k]\). Die neuen Koordinaten im Hauptkomponentenraum sind dann analog \(XA^k=Z^k\), wobei \(Z^k\in \mathbb{R}^{n\times k}\). Insbesondere erhält man für \(k=2\) neu \(n\) Zeilen und zwei Spalten, diese bilden zweidimensionale Koordinaten für unsere \(n\) Objekte im Raum der Hauptkomponenten und können in einer Ebene geplottet werden.
3.2 PCA basierend auf den Worldbank Daten
Wir werden hier nicht versuchen, Variablen zu reduzieren, um schliesslich mit weniger Variablen eine Modellbildung zu machen. Wir stellen uns hier die Frage, ob wir die PCA benutzen können, um mit den im Teil 1 beschriebenen Variablen Gruppen von Ländern zu identifizieren, welche betr. der Variablen bzw. der identifizierten PCAs ähnlich sind.
Wir haben folgende sieben Variablen im Datensatz:
| Name in Original | neuer Name |
|---|---|
| GDP per capita (current US$) | GDP_per_capita |
| Female share of employment in senior and middle management (%) | Females_Management |
| Part time employment, female (% of total female employment) | Female_parttime |
| Proportion of seats held by women in national parliaments (%) | Females_Parliament |
| School enrollment, secondary (% gross) | School_enrollment_secondary |
| Employment to population ratio, 15+, total (%) (national estimate) | Employment |
| Government expenditure on education, total (% of GDP) | Expenditure_education |
head(df_ind_country_nona)## # A tibble: 6 × 8
## Country GDP_per_capita Females_Management Female_parttime Females_Parliament
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Albania 6413. 30.7 26.8 35.7
## 2 Angola 1926. 11.0 34.8 29.5
## 3 Argentina 10738. 37.4 57.5 44.7
## 4 Austria 53649. 34.3 65.0 40.4
## 5 Belarus 7490. 44.0 22.2 40
## 6 Belgium 51656. 34.4 56.2 42
## # ℹ 3 more variables: School_enrollment_secondary <dbl>, Employment <dbl>,
## # Expenditure_education <dbl>Jede Zeile ist entspricht einem Land, was ebenfalls als Spalte aufgeführt ist. Später werden wir noch die Zugehörigkeit zu einer Weltregion hinzufügen.
Explorative Fragestellungen wären dann: Sind sich Länder aus den verschiedenen Weltregionen, z.B. Westeuropa betr. der gewählten Variablen ähnlich?
Wir können nun eine PCA mit unserem Datensatz machen. Wir müssen unsere Daten noch skalieren, bzw. standardisieren, da wir ja im 1. Kapitel gesehen haben, dass die Varianzen sehr unterschiedliche Grössenordnungen haben. Weiter zentrieren wir die Daten um den Nullpunkt.
df_ind_country_nona_num <- as.data.frame(df_ind_country_nona)
rownames(df_ind_country_nona_num) <- df_ind_country_nona_num$Country
df_ind_country_nona_num$Country <- NULL df_pca <- princomp(scale(df_ind_country_nona_num, scale=TRUE, center=TRUE), scores=TRUE)
par(mar=c(6,3,0,0)+.1, las=2)
screeplot(df_pca, main='')
Wir sehen im Plot oben die einzelnen Hauptkomponenten und auf der y-Achse die Varianz (\(\lambda_i\)), welche sie erklären. Bei skalierten Daten gilt \(\sum_i^p \lambda_1 = p\). Hier ist \(p=7\). Die erste Hauptkomponente, welche eine Varianz von ca. \(2.76\) aufweist, erklärt hier also ca. 39 Prozent der Varianz. Hauptkomponente 1 und 2 zusammen erklären ca. 58 Prozent der gesamten Varianz. Besser wäre natürlich mehr. Es ist hier aber auch kein deutlicher Knick im Screeplot vorhanden, welcher uns vorgeben würde, welche Anzahl Hauptkomponenten sinnvoll wäre. Wir beschränken uns auf zwei, da wir das Ziel gesetzt haben, die Daten im zweidimensionalen Raum zu plotten.
eigenwerte <- df_pca$sdev^2
eigenwerte## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## 2.7639491 1.3258682 0.8638077 0.7730536 0.5128386 0.3427883 0.3029404Schauen wir uns noch die Loadings für die ersten zwei Hauptkomponenten an und erstellen einen Plot, wo wir die Loadings gegen die Hauptkomponenten abbilden.
loadings <- df_pca$loadings[,1:2]
loadings <- as.data.frame(loadings)
loadings$Symbol <- row.names(loadings)
loadings <- gather(loadings, 'Component', 'Weight', -Symbol)loadings$Color = loadings$Weight > 0
graph <- ggplot(loadings, aes(x=Symbol, y=Weight, fill=Color)) +
geom_bar(stat='identity', position='identity', width=.75) +
facet_grid(Component ~ ., scales='free_y') +
guides(fill='none') +
ylab('Component Loading') +
theme_bw() +
theme(axis.title.x=element_blank(),
axis.text.x=element_text(angle=90, vjust=0.8))
graph
Wir sehen, dass die erste Hauptkomponente von den Variablen, Expenditure_education, Female_parttime, Females_parliament, GPR_per_capital und school_enrollment_secondary dominiert wird. Es ist eine Art gewichteter Mittelwert dieser Variablen.
Die zweite Hauptkomponente hingegen fängt vor allem die Varianz der Variable Females_Management ein.
Projizieren wir unsere Daten nun einmal auf die ersten beiden Hauptkomponenten.
# Extrahiere die Projektionen der Daten auf die ersten beiden Hauptkomponenten
pc1_pc2_scores <- df_pca$scores[, 1:2]
# Scatterplot der projizierten Daten
plot(pc1_pc2_scores, main="Projektion auf die ersten beiden Hauptkomponenten",
xlab="PC1", ylab="PC2", pch=19, col="blue")
text(pc1_pc2_scores, labels = rownames(pc1_pc2_scores), pos = 3, cex = 0.7, col = "red")
Viele westeuropäische Länder liegen weit links im Plot. Wir untersuchen dies weiter, indem wir unseren Daten noch die Zugehörigkeit zu verschiedenen Weltregionen anfügen. Die Weltregionen wurden gemäss Zuteilung von chatGPT erfasst und sind in einem separaten File gelistet.
regions <- read_csv("liste_country_region2.csv")
pc1_pc2_scores_test <- as.data.frame(pc1_pc2_scores)
pc1_pc2_scores_test$Country <- rownames(pc1_pc2_scores_test)
pca_with_regions_test<-merge(pc1_pc2_scores_test , regions, by = "Country")Nun können wir die Projektion unserer Daten auf die Hauptkomponenten nochmals plotten und die Datenpunkte nach Weltregionen einfärben.
ggplot(pca_with_regions_test, aes(x = Comp.1, y = Comp.2, color=Weltregion)) +
geom_point(size = 4) +
labs(title = "PCA Plot mit Weltregionen",
x = "Principal Component 1",
y = "Principal Component 2") +
scale_color_brewer(palette = "Set1") +
theme_minimal() 
Wir sehen, dass sie die westeuropäischen Länder auf der linken Seite von 0 der ersten Hauptkomponente bewegen. Es sind also diejenigen Länder, die bei den Variablen Expenditure_education, Female_parttime, Females_parliament, GPR_per_capita und school_enrollment_secondary, welche die erste Hauptkomponenten massgeblich (positiv) beeinflussen, eher höhere Werte ausweisen. Im Bereich der zweiten Hauptkomponenten bewegen sie sich alle eher im mittleren Band (der Werte aller Länder). Die meisten zwischen 1 und -1. D.h. betr. der Variable Females_management, welche die zweite Hauptkomponenten dominierte, bewegen sie sich eher im Mittelfeld.
Die Länder der Region Mittel- und Südosteuropa clustern demgegenüber eher auf der rechten Seite von 0 der PC1 zusammen. Sie weisen vermutlich bei den massgebenden Variablen eher mittlere bis tiefere Werte auf. Betreffend der zweiten Hauptkomponente liegen sie dichter als die Länder von Westeuropa zusammen, d.h. im Bereich Females_management liegen sie vermutlich insgesamt dicht beieinander und bewegen sich im Mittelfeld.
Die weiteren Länder clustern weniger schön nach Regionen zusammen. Für die beiden Länder aus Subsahara-Afrika kann man jedoch sagen, dass sie sich nur betr. der zweiten Hauptkomponente unterscheiden. Die Länder aus Asien und Pazifik liegen betr. der ersten Hauptkomponenten ebenfalls relativ nahe beieinander aber weisen betr. der zweiten Hauptkomponente eine grosse Varianz auf. Die Länder aus Lateinamerika hingegen liegen betr. der zweiten Hauptkomponente nahe beieinander und weisen dafür betr. der ersten Hauptkomponente eine grosse Varianz auf. Wenn man wieder die Anteile der verschiedenen Variablen an den Hauptkomponenten betrachtet, dann heisst das, dass die Länder aus Lateinamerika vermutlich ähnliche Werte bei Females_Management aufweisen, aber betr. der dominierenden Variablen von PC1 grössere Differenzen aufweisen.
Die Zugehörigkeit zu verschiedenen Weltregionen hat offenbar einen Zusammenhang mit den gemessenen Variablen/Indikatoren der einzelnen Länder. Die Farben der Regionen sind auf dem zweidimensionalen Plot nicht völlig zufällig verteilt, sondern verschiedene Länder einer Region liegen eher näher zusammen und man kann verschiedene Gruppen identifizieren. Insbesondere sieht man, dass PC1 zwischen der Region “Mittel- und Südosteuropa” und der Region “Westeuropa” differenziert.
Man kann also daraus schon einmal die Hypothese ableiten, dass Länder aus Westeuropa betr. der Variablen Expenditure_education, Female_parttime, Females_parliament, GPR_per_capita und school_enrollment_secondary höhere Werte haben als Länder aus Mittel- und Südosteuropa. Dies könnte man dann mit weiteren statistischen Tests verifizieren.
Hier wollen wir aber weiter untersuchen, ob eine andere Methode unsere Daten vielleicht ebenfalls oder sogar noch besser nach Weltregionen aufteilen kann.
4 Clusteranalyse
4.1 Clusteranalyse - Theoretische Grundlagen
(Der Theorieteil basiert auf (Handl and Kuhlenkasper 2017) und den Schulungsunterlagen, i.e. YouTube-Videos von DATAtab und Prof. Victor Lavrenko (Univ. Edinburgh).)
Bei der Clusteranalyse wird versucht, mit geeigneten Methoden ähnliche Datenpunkte/Objekte in einem Datensatz zu aggregieren, i.e. zu Clustern zusammenzufügen, und so unterliegende Strukturen und Hierarchien in den Daten aufzudecken.
Es gibt zwei Hauptansätze wie man Datensätze clustern kann, einerseits die hierarchischen Clustering-Verfahren und andererseits die partitionierenden Clustering-Verfahren. Die hierarchischen Verfahren kann man in agglomerative und divisive Verfahren unterteilen. Bei ersteren bildet mit einem Bottom-up Approach zunächst jeder Punkt einen eigenen Cluster, danach werden Cluster immer weiter zusammengefügt bis ein vordefiniertes Stopp-Kriterium erreicht ist oder bis alle Punkte in einem Cluster landen. Bei letzteren werden aus einem einzigen anfänglichen Cluster, welcher alle Datenpunkte enthält, durch Unterteilung immer weitere Verfeinerungen geschafft (Top-down Approach).
Bei den partitionierenden Verfahren legt man die Anzahl Cluster im Vornherein fest und startet mit zufällig gewählten Clustern. Danach verbessert man diese mit einem Algorithmus bis ebenfalls ein Stopp-Kriterium erreicht wird.
Zunächst müssen wir aber definieren, was wir unter Abständen zwischen zwei Datenpunkten/Objekten und schliesslich verallgemeinert zwischen zwei Clustern verstehen. Dafür brauchen wir Metriken.
4.1.1 Metriken
Eine Metrik gibt uns eine Möglichkeit, Abstände zwischen zwei Punkten zu messen.
Sei \(X\) eine beliebige Menge. Eine Funktion
\[ d: X\times X \rightarrow \mathbb{R}_0^{+} \] heisst Metrik auf \(X\), wenn für beliebige \(x,y,z \in X\) folgendes gilt:
- \(d(x,y)=d(y,x)\)
- \(d(x,y)=0 \Leftrightarrow x=y\)
- \(d(x,y) \leq d(x,z)+d(z,y)\)
(vgl. https://de.wikipedia.org/wiki/Metrischer_Raum).
Schauen wir uns ein paar Beispiele für Metriken auf \(\mathbb{R}^n\) an, die im Folgenden auch benutzt werden. Seien die Punkte \(x=(x_1, x_2,...,x_n)^T\) und \(y=(y_1,y_2, ..., y_n)^T\) aus \(\mathbb{R}^n\):
Euklidische Distanz: \[ d(x,y)=\left(\sum_{i=1}^{n}(x_i - y_i)^2\right)^\frac{1}{2} \]
Manhattan Distanz: \[d(x,y)=\sum_{i=1}^n|x_i-y_i|\]
Maximum (Tschebyschew) Distanz : \[d(x,y)=max_{i}|x_i-y_i|\]
Mahalanobis Distanz: \[d(x,y)=(x-y)^T \Sigma^{-1}(x-y)\text{, mit der Kovarianz-Matrix }\Sigma\]
4.1.2 Linkage-Verfahren
Nun müssen wir noch Distanzen zwischen Clustern definieren. Sei also ein Cluster \(C_i\) eine Menge von Punkten \(c_{ki}\) für \(k=1,..,p_i\). Die Distanz zwischen zwei Clustern \(C_i\) und \(C_j\) kann mit folgenden Methoden definiert werden:
- Single Linkage: \[ \min_{c_{ki}\in C_i, c_{lj}\in C_j} d(c_{ki},c_{lj}) \]
- Complete Linkage: \[ \max_{c_{ki}\in C_i, c_{lj}\in C_j} d(c_{ki},c_{lj}) \]
- Average Linkage: \[ \frac{1}{|C_i|\cdot |C_j|}\sum_{c_{lj}\in C_j} \sum_{c_{ki}\in C_i} d(c_{ki},c_{lj}) \]
- Ward Linkage: \[ \frac{|C_i|\cdot |C_j|}{|C_i|+|C_j|} d_{Euklid}(\bar{x}_{C_i}, \bar{x}_{C_j})^2 \]
\[ \text{wobei }\bar{x}_{C_i},\bar{x}_{C_j}\text{ die Schwerpunkte (Mittelwerte) der Cluster sind.} \]
Wir haben nun also Methoden, um Abstände zwischen Punkten, sowie zwischen Clustern zu berechnen. Ein Abstand eines einzelnen Punkts zu einem Cluster kann ebenfalls berechnet werden, indem der Punkt als Cluster mit nur einem Objekt gesehen wird.
4.1.3 Hierarchisches Clustering
Wir schauen uns hier das hierarchische agglomerative Clustering als Pseudoalgorithmus an:
- Bestimme eine zu verwendende Metrik und das Verfahren zur Bestimmung der Distanz zweier Cluster
- Jedes Objekt unserer Daten bildet einen eigenen Cluster
- Finde in der Menge aller Cluster die zwei Cluster mit der minimalen Distanz und füge die beiden Cluster zu einem einzigen Cluster zusammen.
- Wiederhole Schritt 3 bis zum Erreichen eines Stopp-Kriteriums (z.B. vordefinierte Distanz / Anzahl Cluster / alle Punkte in einem Cluster)
Man kann bei dieser Vorgehensweise ein sogenanntes Dendrogramm (Baumdiagramm) erstellen, welches den Verlauf des Algorithmus visualisert und uns Einblick in die Struktur der Daten gibt. Auf der x-Achse liegen alle Objekte auf der y-Achse wird jeweils der Abstand aufgetragen. Wenn zwei Objekte oder Cluster zu einem neuen Cluster zusammengefügt werden, wird deren Distanz in senkrechter Richtung abgetragen. Wenn man das Dendrogramm auf einer beliebigen Höhe \(y_\text{cutoff}\) durchschneidet, erhält man die Cluster zu dieser Distanz, d.h. alle resultierenden Cluster haben zueinander eine Distanz, welche mindestens \(y_\text{cutoff}\) beträgt oder grösser ist. Beispiele zum Dendrogramm, bzw. die grafische Darstellung davon, werden wir direkt anhand unserer Daten anschauen.
4.1.4 Partitionierendes Clustering: \(k\)-means
Bei \(k\)-means handelt es sich um ein sogenanntes partitionierendes Clustering-Verfahren. Die Anzahl Cluster wird im Vornherein festgelegt. Wir bestimmen \(k\) als Anzahl der Cluster, \(n\) ist wieder die Anzahl Datenpunkte, bzw. Objekte, in unserem Datensatz. Der \(k\)-means-Algorithmus funktioniert folgendermassen:
- \(k\) Datenpunkte/Objekte aus dem Datensatz werden zufällig ausgewählt und als Clusterzentren festgelegt.
- Alle Objekte werden zu dem Cluster zugeordnet, zu welchem Clusterzentrum sie am nächsten sind.
- Die Clusterzentren werden neu als Schwerpunkte/Mittelwerte der Punkte/Objekte in einem Cluster festgelegt.
- Schritt 2 und 3 werden so lange wiederholt, bis sich nichts mehr verändert oder ein Stopp-Kriterium erfüllt ist.
Es ist offensichtlich, dass diese Methode von den initial gewählten Clustern abhängt oder zumindest abhängen kann. Deswegen wird dieses Verfahren üblicherweise mehrfach mit verschiedenen Startpunkten durchgeführt und am Schluss das “beste” Ergebnis gewählt. Die Güte eines Ergebnisses wird dabei durch die Summe aller Abstände aller Punkte zu ihren jeweiligen Clusterzentren gemessen, diese soll möglichst minimal sein.
Ein zweites Problem ist die Wahl der Anzahl Cluster. Wie können wir überprüfen, ob unsere Wahl für \(k\) sinnvoll war oder ob wir allenfalls für unsere Daten eine andere, bessere Wahl treffen müssen. Dazu wird normalerweise die Ellbogen-Methoden verwendet. Hierbei wird \(k\) auf der \(x\)-Achse aufgetragen und die Summe der quadrierten Abstände aller Punkte zu ihren jeweiligen Clusterzentren auf der \(y\)-Achse. Die Summe der quadrierten Abstände nimmt mit wachsendem \(k\) stetig ab, aber diese Kurve macht in der Regel bei einem gewissen \(k\) einen Knick. Dieses \(k\) wird dann als optimal gewählt.
4.2 Clusteranalyse basierend auf den Worldbank Daten
Wir fahren mit den gleichen Daten fort wie bei der PCA und beschäftigen uns mit den gleichen Fragestellungen.
Wir wollen also auch hier herausfinden, ob sich die Länder aus den verschiedenen Weltregionen, z.B. Westeuropa, betr. der gewählten Variablen - Ausgaben für die Bildung, “School Enrollment”, Beschäftigungsrate, Anteil weiblicher Mitglieder im Parlament, Anteil weiblicher Angestellter im höheren und mittleren Management, Anteil weiblicher Teilzeitarbeit und BIP pro Kopf - ähnlich sind. Wenn wir die Daten nun clustern, landen ähnliche Länder in gleichen Clustern. Wir können uns dann ansehen, ob die Cluster die Weltregionen separieren.
Nehmen wir also wieder unseren Datensatz mit den 61 Ländern und fügen die Spalte “Weltregion” hinzu.
df_ind_country_nona_region <- merge(df_ind_country_nona, regions, by = "Country")head(df_ind_country_nona_region)## Country GDP_per_capita Females_Management Female_parttime
## 1 Albania 6413.283 30.675 26.85
## 2 Angola 1925.875 11.014 34.78
## 3 Argentina 10738.018 37.360 57.53
## 4 Austria 53648.719 34.309 65.03
## 5 Belarus 7489.719 43.983 22.16
## 6 Belgium 51655.788 34.446 56.22
## Females_Parliament School_enrollment_secondary Employment
## 1 35.71429 97.66503 52.836
## 2 29.54545 52.54377 64.306
## 3 44.74708 116.55945 55.085
## 4 40.43716 102.30766 57.202
## 5 40.00000 95.21339 67.307
## 6 42.00000 143.16351 51.079
## Expenditure_education Weltregion
## 1 3.022560 Mittel- und Südosteuropa
## 2 2.297109 Subsahara-Afrika
## 3 4.641170 Lateinamerika
## 4 5.494110 Westeuropa
## 5 4.606410 Osteuropa und Zentralasien
## 6 6.355520 Westeuropa
## World_Region
## 1 Middle and Southern Europe
## 2 Sub-Saharan Africa
## 3 Latin America and the Carribean
## 4 Western Europe
## 5 Eastern Europe and Central Asia
## 6 Western EuropeDamit die Cluster nicht von einer Variable und deren Varianz dominiert werden, müssen wir wiederum skalieren.
df_ind.scaled<-cbind(scale(df_ind_country_nona_region[, 2:8]), Weltregion=df_ind_country_nona_region$Weltregion, Land=df_ind_country_nona_region$Country,Land_Reg=df_ind_country_nona_region$Land_Reg)4.2.1 Hierarchisches Clustering
Wir berechnen die Distanzmatrix mit der euklidischen Metrik und clustern danach hierarchisch mit dem Verfahren “Average Linkage”.
dist.df_ind_scaled <- dist(df_ind.scaled[,1:7], method ="euclidean")
fit.ave <-hclust(dist.df_ind_scaled, method = "average")
plot(fit.ave, hang=-1, labels=df_ind_country_nona_region$Country, cex=.7,main = "Dendrogramm - Metrik: Euklidisch, Linkage: Average", xlab="")
Wir haben acht verschiedene Weltregionen und können nun den Cut-off beim Clustering einmal so setzen, dass wir acht Cluster erhalten. D.h. wir schneiden das obige Dendrogramm genau bei der Höhe durch, dass darunter exakt acht Cluster übrig bleiben.
Diese jeweiligen Cluster schreiben wir noch in unseren Datensatz und plotten dann wieder das Dendrogramm mit den markierten Clustern. Dabei visualisieren wir auch gerade die Aufteilung nach Weltregionen.
# Zuordnung von 8 Clustern zu unserem originalen Datensatz
df_ind_country_nona_region$clus.ave <- cutree(fit.ave, k=8)# 1. Dendrogramm erzeugen
dend <- as.dendrogram(fit.ave)
# 2. Länder in Dendrogramm-Reihenfolge
leaf_order <- order.dendrogram(dend)
country_names <- df_ind_country_nona_region$Country[leaf_order]
# 3. Weltregionen passend zur Dendrogrammreihenfolge
weltregionen_ordered <- df_ind_country_nona_region$Weltregion[leaf_order]
weltregionen_ordered <- factor(weltregionen_ordered)
# 4. Farben pro Region
farben <- brewer.pal(n = length(levels(weltregionen_ordered)), "Set1")
label_colors <- farben[weltregionen_ordered]
# 5. Länderlabels und Farben setzen
labels(dend) <- country_names
labels_colors(dend) <- label_colors
# 6. Dendrogramm plotten
op <- par(mar = c(10, 2, 2, 2)) # mehr Platz rechts
plot(dend, cex = 0.6, main = "Dendrogramm - Metrik: Euklidisch, Linkage: Average")
rect.dendrogram(dend, k = 8, border = "red")
legend("topright",
legend = levels(weltregionen_ordered),
fill = farben,
border = NA,
bty = "n",
cex = 0.7,
ncol = 2,
inset = c(-0.02, -0.05),
xpd = TRUE) 
par(op)Wir können uns noch in Tabellenform anschauen, wie viele Länder der verschiedenen Weltregionen in den einzelnen Clustern enden.
table(df_ind_country_nona_region$Weltregion,df_ind_country_nona_region$clus.ave)##
## 1 2 3 4 5 6 7 8
## Asien und Pazifik 8 2 0 0 0 0 0 0
## Lateinamerika 3 0 3 1 0 0 0 0
## Mittel- und Südosteuropa 17 0 0 0 0 0 0 0
## Nordafrika und Nahost 0 0 1 0 0 1 0 0
## Nordamerika 0 0 1 0 0 0 0 0
## Osteuropa und Zentralasien 4 0 0 0 0 0 0 0
## Subsahara-Afrika 0 1 0 0 0 0 1 1
## Westeuropa 1 0 14 0 2 0 0 0Wir sehen, dass mit dieser Clustering-Methoden die westeuropäischen Länder alle ausser Luxemburg, Irland und Italien dem Cluster 3 zugeordnet werden. Diesem Cluster sind weiter nur wenige andere Länder zugeordnet. Die Region Westeuropa wird durch das Clustering also relativ gut vom Rest der Daten separiert. Jedoch sind die weiteren Cluster weniger aussagekräftig. Es gibt noch einen weiteren sehr grossen Cluster (Cluster 1), welcher betr. Weltregionen durchmischt ist, aber alle Ländern aus Mittel- und Südosteuropa enthält und von diesen dominiert wird. Daneben gibt es sechs weitere Cluster, welche alle höchstens 3 Länder beinhalten.
Versuchen wir es also noch mit einer anderen Methode. Wir nehmen als Metrik wieder die euklidische Metrik, für die Distanzberechnung zwischen Clustern nehmen wir jedoch das Ward-Verfahren.
fit.ward <-hclust(dist.df_ind_scaled, method = "ward.D2")
plot(fit.ward, hang=-1, labels=df_ind_country_nona_region$Weltregion, cex=.7,main = "Dendrogramm - Metrik: Euklidisch, Linkage: Ward")
# 8-Cluster
df_ind_country_nona_region$clus.ward <- cutree(fit.ward, k=8)# 1. Dendrogramm erzeugen
dend <- as.dendrogram(fit.ward)
# 2. Länder in Dendrogramm-Reihenfolge
leaf_order <- order.dendrogram(dend)
country_names <- df_ind_country_nona_region$Country[leaf_order]
# 3. Weltregionen passend zur Dendrogrammreihenfolge
weltregionen_ordered <- df_ind_country_nona_region$Weltregion[leaf_order]
weltregionen_ordered <- factor(weltregionen_ordered)
# 4. Farben pro Region
farben <- brewer.pal(n = length(levels(weltregionen_ordered)), "Set1")
label_colors <- farben[weltregionen_ordered]
# 5. Länderlabels und Farben setzen
labels(dend) <- country_names
labels_colors(dend) <- label_colors
# 6. Dendrogramm plotten
op <- par(mar = c(10, 2, 2, 2))
plot(dend, cex = 0.6, main = "Dendrogramm - Metrik: Euklidisch, Linkage: Ward")
rect.dendrogram(dend, k = 8, border = "red")
legend("topright",
legend = levels(weltregionen_ordered),
fill = farben,
border = NA,
bty = "n",
cex = 0.7,
ncol = 2,
inset = c(-0.04, 0.05),
xpd = TRUE) 
par(op)table(df_ind_country_nona_region$Weltregion,df_ind_country_nona_region$clus.ward)##
## 1 2 3 4 5 6 7 8
## Asien und Pazifik 2 5 0 0 2 1 0 0
## Lateinamerika 2 0 3 1 0 1 0 0
## Mittel- und Südosteuropa 4 0 0 0 4 9 0 0
## Nordafrika und Nahost 1 0 0 0 0 1 0 0
## Nordamerika 0 0 0 0 0 1 0 0
## Osteuropa und Zentralasien 3 0 0 0 0 1 0 0
## Subsahara-Afrika 0 1 1 1 0 0 0 0
## Westeuropa 0 0 8 0 1 0 6 2Auch hier sehen wir, dass die westeuropäischen Länder sich von den restlichen Ländern separieren. Die Cluster 7 und 8 bestehen nur aus westeuropäischen Ländern, Cluster 3 wird von westeuropäischen Ländern dominiert. Auf der einer höheren Distanz-Ebene würden diese Cluster zusammenschmelzen. Sie sind in diesem Sinne näher zueinander und weiter entfernt von den anderen Clustern. Weiter sehen wir, dass die Länder aus Mittel- und Südosteuropa den Clustern 1, 5 und 6 zugeordnet werden und diese hauptsächlich ausmachen. Cluster 2 besteht hauptsächlich aus Ländern aus der Region Asien und Pazifik.
4.2.2 \(k\)-means Clustering
Wir können auch noch versuchen, mit dem \(k\)-means Algorithmus zu clustern. Standardmässig wird hier die euklidische Distanz für die Berechnung des Abstands zu den Clusterzentren verwendet.
Wir wählen \(k=8\), also acht Cluster, weil wir acht Weltregionen haben.
Wir können uns aber noch die Fehlerquadratsumme gegen \(k\), die Anzahl Cluster, anschauen.
# Berechne Fehlerquadratsumme
wss.fitmeans <- (nrow(df_ind.scaled[,1:7])-1)*sum(apply(df_ind.scaled[,1:7],2,var))
for( i in 2:15){
wss.fitmeans[i] <- sum(kmeans(df_ind.scaled[,1:7], centers=i)$withinss)
}
# Fehlerquadratsumme gegen Anzahl Cluster plotten
plot(1:15, wss.fitmeans, type="b", xlab="Anzahl Cluster",
ylab="Fehlerquadratsumme")
Wir können aus dem Plot Fehlerquadratsumme gegen Anzahl Cluster keinen klaren Knick bei einem bestimmten Wert von \(k\) erkennen. In diesem Sinne spricht also nichts gegen acht Cluster.
# 8 Cluster-Lösung mit 8 zufälligen initialen Clusterzentren
fit.means <- kmeans(df_ind.scaled[,1:7], 8, nstart=25)
# Clusterzugehörigkeit speichern
df_ind_country_nona_region$clus.means <- fit.means$clusterWir können auch hier eine Art Dendrogramm plotten. Dabei wird jeweils ein Baumdiagramm mit Hilfe der Distanzen zwischen den Clusterzentren erstellt.
# 1. Clusterzentren auswählen und Dendrogramm berechnen
centers <- fit.means$centers
hc_centers <- hclust(dist(centers))
# 2. Cluster-Zuordnung & Länderinfos
cluster_assignments <- fit.means$cluster
country_names <- df_ind_country_nona_region$Country
weltregionen <- df_ind_country_nona_region$Weltregion
# 3. Farben gemäss Weltregionen
weltregionen_factor <- factor(weltregionen)
region_colors <- setNames(brewer.pal(length(levels(weltregionen_factor)), "Set1"),
levels(weltregionen_factor))
point_colors <- region_colors[weltregionen]
# 4. Gruppiere die Länder nach Cluster
cluster_country_indices <- lapply(1:max(cluster_assignments), function(k) {
which(cluster_assignments == k)
})
# 5. Margins setzen
old_par <- par(no.readonly = TRUE)
par(mar = c(6, 4, 4, 8) + 0.1, xpd = TRUE)
# 6. Dendrogramm plotten mit angepasstem ylim
max_height <- max(hc_centers$height)
plot(hc_centers,
labels = FALSE,
main = "Baumdiagramm der Clusterzentren von k-means",
xlab = "",
ylim = c(-0.1 * max_height, max_height)) # y-Achse erweitert
# 7. Punkte unterhalb der Cluster zeichnen
x_positions <- 1:length(cluster_country_indices)
for (i in seq_along(x_positions)) {
idxs <- cluster_country_indices[[hc_centers$order[i]]]
y_val <- -0.02 * max_height # Punkte näher am Dendrogramm
x_jitter <- seq(-0.2, 0.2, length.out = length(idxs))
points(x = rep(x_positions[i], length(idxs)) + x_jitter,
y = rep(y_val, length(idxs)),
col = point_colors[idxs],
pch = 16, cex = 1.7)
}
# 8. Legende
legend("topright",
legend = names(region_colors),
col = region_colors,
pch = 16,
pt.cex = 1.5,
border = NA,
bty = "n",
cex = 0.8,
inset = c(-0.2, 0))
# 9. Grafikparameter zurücksetzen
par(old_par)Die einzelnen Länder werden, als nach Weltregion eingefärbte Punkte, unter den jeweiligen Clustern dargestellt.
table( df_ind_country_nona_region$Weltregion, df_ind_country_nona_region$clus.means)##
## 1 2 3 4 5 6 7 8
## Asien und Pazifik 0 0 2 2 3 0 0 3
## Lateinamerika 1 0 2 0 0 3 0 1
## Mittel- und Südosteuropa 0 0 4 5 0 0 0 8
## Nordafrika und Nahost 0 0 1 0 0 0 0 1
## Nordamerika 0 0 0 0 0 0 0 1
## Osteuropa und Zentralasien 0 0 3 0 0 0 0 1
## Subsahara-Afrika 1 0 0 0 1 1 0 0
## Westeuropa 0 10 0 1 0 3 3 0Wir sehen, dass mit dieser Methode zwei Cluster entstehen, welche nur westeuropäische Länder enthalten. Weitere vier westeuropäische Länder sind in zwei gemischten Cluster vertreten. Die mittel- und südosteuropäischen Länder sind in drei verschiedenen Clustern vertreten. Man sieht also auch hier, dass die Mehrheit der westeuropäischen Länder von den restlichen Ländern durch das Clustering separiert wird.
5 Fazit Clustering und PCA
Die Resultate der Clusterings sowohl hierarchisch als auch mit \(k\)-means entsprechen dem Resultat der PCA, wo wir bereits gesehen haben, dass in der Ebene von PCA 1 und PCA 2 sich die Länder aus Westeuropa und die Länder aus Mittel- und Südosteuropa relativ gut gruppieren und sich insbesondere die Länder aus Westeuropa von den restlichen Ländern separieren. Dasselbe sehen wir auch hier. Westeuropa separiert sich bei allen angewendeten Verfahren relativ gut von den Ländern aus den restlichen Regionen. Die Länder aus Mittel- und Südosteuropa gruppieren sich ebenfalls, sind jedoch nicht so klar abgrenzbar.
Wir können also schliessen, dass die Länder aus Westeuropa und teilweise bzw. in schwächerer Form auch die Länder aus Mittel- und Südosteuropa sich bezüglich der gewählten Variablen ähnlich sind. Wir können dank der PCA auch Aussagen darüber treffen, bei welchen Variablen sie sich ähnlich sind und ob sie dort eher extremere Werte oder Werte im Mittelfeld annehmen. Dazu verweise ich auf meine Ausführungen im Kapitel zu PCA. Nur mit dem Clustering hätte man solche Aussagen nicht treffen können, es ist also sinnvoll die verschiedenen Methoden miteinander zu kombinieren.
5.1 Überprüfung mit MANOVA
Ob sich die Weltregionen betr. der hier untersuchten Variablen wirklich signifikant unterscheiden, kann man nun mit MANOVA weiter untersuchen.
# 1. Die relevanten Spaltennamen definieren (abhängige Variablen)
vars <- c("GDP_per_capita", "Females_Management", "Female_parttime",
"Females_Parliament", "School_enrollment_secondary",
"Employment", "Expenditure_education")
# 2. Nur die relevanten Spalten herausholen und skalieren
scaled_vars <- scale(df_ind_country_nona_region[, vars])
# 3. Gruppierungsvariable (Weltregion) extrahieren
group <- as.factor(df_ind_country_nona_region$Weltregion)
# 4. MANOVA-Modell berechnen
manova_model <- manova(scaled_vars ~ group)
# 5. Zusammenfassung der multivariaten Tests (z. B. Wilks' Lambda)
summary(manova_model, test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## group 7 0.037389 4.5156 49 243.03 1.7e-15 ***
## Residuals 53
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ein kleines Wilks’ Lambda (hier 0.037389) bedeutet, dass die Gruppenzugehörigkeit (Weltregion) einen grossen Anteil der Varianz im multivariaten Raum erklärt. Zudem ist der P-Wert für die approximierte F-Statistik sehr klein. Das heisst, die Weltregionen unterscheiden sich betr. dieser Variablen zusammengenommen signifikant voneinander.
Die MANOVA-Analyse stützt also die Schlussfolgerungen, welche wir aus der PCA und der Cluster-Analyse gezogen haben.
Hilfsmittelverzeichnis
| Welches Hilfsmittel wurde eingesetzt? | Wozu wurde es eingesetzt? | Betroffene Stellen |
|---|---|---|
| ChatGPT | zwecks Findung von Methoden in R | einzelne Codestellen in der gesamten Arbeit (e.g. slice()) |
| ChatGPT | Erstellung eines Files mit Zuordnung von Ländern zu Weltregionen | PCA |
| Kursunterlagen | Theorie betr. Mult. lin Reg., PCA und Clusteranalyse | Gesamte Arbeit |
| Kursunterlagen, Tutorial/Fallstudien aus diesem Kurs | Methoden in R für Mult. lin Reg., PCA und Clusteranalyse | Mult. lin. Reg., PCA und Clusteranalyse |
| Kursunterlagen und meine Semesterarbeiten aus CAS Grundlagen Data Science | Methoden in R | e.g. Verwendung von grid.arrange |
| ChatGPT | Erstellen von Plots | Code mit Kommentaren zum Plotten von Punkten in der PCA-Ebene (PCA) |
| ChatGPT | Erstellen von Plots | Code mit Kommentaren zum Plotten von Dendrogrammen mit nach Weltregionen gefärbten Labels (Clusteranalyse) |
| ChatGPT | Code mit Kommentaren für Manova | Fazit PCA und Clusteranalyse - MANOVA |
| ChatGPT | Erklärung wie der Output von MANOVA interpretiert werden muss | Fazit PCA und Clusteranalyse - MANOVA |
| ChatGPT | Kontrolle Rechtschreibung (ohne Stil) | Gesamtes File |
Eigenständigkeitserklärung
Mit der Abgabe dieser Arbeit bestätige ich, dass ich die vorliegende Arbeit selbstständig verfasst habe, dass alle sinngemäss und wörtlich übernommenen Textstellen aus fremden Quellen kenntlich gemacht wurden, dass alle mit Hilfsmitteln erbrachten Teile der Arbeit präzise deklariert wurden, dass keine anderen als die im Hilfsmittelverzeichnis aufgeführten Hilfsmittel verwendet wurden, dass das Thema, die Arbeit oder Teile davon nicht bereits Gegenstand eines Leistungsnachweises eines anderen Moduls waren, sofern dies nicht ausdrücklich mit der Referentin oder dem Referenten im Voraus vereinbart wurde, dass ich mir bewusst bin, dass meine Arbeit elektronisch auf Plagiate und auf Drittautorschaft menschlichen oder technischen Ursprungs überprüft werden kann und ich hiermit der FFHS das Nutzungsrecht so weit einräume, wie es für diese Verwaltungshandlungen notwendig ist.
Küsnacht, 15.06.2025
